Load the required libraries.
library(tidyverse)
library(sf)
library(here)
library(readxl)
library(scales)
library(DT)
library(brms)
library(tidybayes)
library(patchwork)
library(marginaleffects)
library(ggrepel)
library(scico)
library(ggdensity)
library(ggpubr)
library(ggsn)
Functions that we will use throughout the script
#labeller for years
year_labels <- c(1950:1963)
#The Glasgow mass minuture chest X-ray campaign happened between 11th March and 12th April 1957
#Segment for graphs to match ACF period
acf_start <- decimal_date(ymd("1957-03-11"))
acf_end <- decimal_date(ymd("1957-04-12"))
Function for counterfactual plots
plot_counterfactual <- function(model_data, model, population_denominator, outcome, grouping_var=NULL, ...){
#labeller for years
year_labels <- c(1950:1963)
#The Glasgow mass minuture chest X-ray campaign happened between 11th March and 12th April 1957
#Segment for graphs to match ACF period
acf_start <- decimal_date(ymd("1957-03-11"))
acf_end <- decimal_date(ymd("1957-04-12"))
summary <- {{model_data}} %>%
select(year, year2, y_num, acf_period, {{population_denominator}}, {{outcome}}, {{grouping_var}}) %>%
add_epred_draws({{model}}) %>%
group_by(year2, acf_period, {{grouping_var}}) %>%
median_qi() %>%
mutate(.epred_inc = .epred/{{population_denominator}}*100000,
.epred_inc.lower = .epred.lower/{{population_denominator}}*100000,
.epred_inc.upper = .epred.upper/{{population_denominator}}*100000) %>%
mutate(acf_period = case_when(acf_period=="a. pre-acf" ~ "Before Intervention",
acf_period=="c. post-acf" ~ "Post Intervention"))
#create the counterfactual (no intervention), and summarise
counterfact <-
add_epred_draws(object = {{model}},
newdata = {{model_data}} %>%
select(year, year2, y_num, {{population_denominator}}, {{grouping_var}}, {{outcome}}) %>%
mutate(acf_period = "a. pre-acf")) %>%
group_by(year2, acf_period, {{grouping_var}}) %>%
median_qi() %>%
mutate(.epred_inc = .epred/{{population_denominator}}*100000,
.epred_inc.lower = .epred.lower/{{population_denominator}}*100000,
.epred_inc.upper = .epred.upper/{{population_denominator}}*100000) %>%
mutate(acf_period = case_when(acf_period=="a. pre-acf" ~ "Before Intervention",
acf_period=="c. post-acf" ~ "Post Intervention"))
#plot the intervention effect
p <- summary %>%
droplevels() %>%
ggplot() +
geom_ribbon(aes(ymin=.epred_inc.lower, ymax=.epred_inc.upper, x=year2, group = acf_period, fill=acf_period), alpha=0.5) +
geom_ribbon(data = counterfact %>% filter(year>=1956),
aes(ymin=.epred_inc.lower, ymax=.epred_inc.upper, x=year2, fill="Counterfactual"), alpha=0.5) +
geom_line(data = counterfact %>% filter(year>=1956),
aes(y=.epred_inc, x=year2, colour="Counterfactual")) +
geom_line(aes(y=.epred_inc, x=year2, group=acf_period, colour=acf_period)) +
geom_point(data = {{model_data}}, aes(y={{outcome}}, x=year2, shape=acf_period), size=2) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
theme_ggdist() +
scale_y_continuous(labels=comma) +
scale_x_continuous(labels = year_labels,
breaks = year_labels) +
scale_fill_manual(values = c("#DE0D92", "grey50", "#4D6CFA") , name="") +
scale_colour_manual(values = c("#DE0D92", "grey50", "#4D6CFA") , name="") +
scale_shape_discrete(name="") +
labs(
x = "Year",
y = "Case notification rate (per 100,000)",
caption = "Mass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)"
) +
theme(legend.position = "bottom",
panel.border = element_rect(colour = "grey78", fill=NA),
title = element_text(size=14),
axis.text.x = element_text(size=10, angle = 90, hjust=1, vjust=0.5),
legend.text = element_text(size=12)) +
guides(shape="none")
facet_vars <- vars(...)
if (length(facet_vars) != 0) {
p <- p + facet_wrap(facet_vars)
}
p
}
Function for calculating measures of change over time
summarise_change <- function(model_data, model, population_denominator, grouping_var=NULL){
#a. immediate change
nd_immediate <- {{model_data}} %>%
filter(year %in% c(1956:1957)) %>%
select(acf_period, year, y_num, {{population_denominator}}, {{grouping_var}})
#Calcuate incidence per draw, then summarise.
immediate_change <- add_epred_draws({{model}},
newdata=nd_immediate) %>%
mutate(epred_inc100k = .epred/{{population_denominator}}) %>%
group_by(.draw, {{grouping_var}}) %>%
mutate(acf_inc100k_diff = last(epred_inc100k)-first(epred_inc100k),
acf_inc100k_rr = last(epred_inc100k)/first(epred_inc100k)) %>%
ungroup() %>%
group_by({{grouping_var}}) %>%
mean_qi(acf_inc100k_diff, acf_inc100k_rr) %>%
mutate(change = "Immediate change") %>%
ungroup()
#b. post-ACF change
nd_post <- {{model_data}} %>%
filter(year %in% c(1956,1958)) %>%
select(acf_period, year, y_num, {{population_denominator}}, {{grouping_var}})
#Calcuate incidence per draw, then summarise.
post_change <- add_epred_draws({{model}},
newdata=nd_post) %>%
mutate(epred_inc100k = .epred/{{population_denominator}}) %>%
group_by(.draw, {{grouping_var}}) %>%
mutate(acf_inc100k_diff = last(epred_inc100k)-first(epred_inc100k),
acf_inc100k_rr = last(epred_inc100k)/first(epred_inc100k)) %>%
ungroup() %>%
group_by({{grouping_var}}) %>%
mean_qi(acf_inc100k_diff, acf_inc100k_rr) %>%
mutate(change = "Post-ACF change") %>%
ungroup()
#c. change in slope post vs. pre-ACF
slope_change <- {{model_data}} %>%
select(year, year2, y_num, acf_period, {{population_denominator}}, {{grouping_var}}) %>%
filter(year!=1957) %>%
add_epred_draws({{model}}) %>%
mutate(inc_100k = .epred/{{population_denominator}}*100000) %>%
group_by(year, {{grouping_var}}, acf_period, ) %>%
mean_qi(inc_100k) %>%
ungroup() %>%
mutate(n_years = length(year), .by=c(acf_period, {{grouping_var}})) %>%
summarise(pct_change_epred_overall = (((last(inc_100k) - first(inc_100k))/first(inc_100k))),
pct_change_lower_overall = (((last(.lower) - first(.lower))/first(.lower))),
pct_change_upper_overall = (((last(.upper) - first(.upper))/first(.upper))),
pct_change_epred_annual = (((last(inc_100k) - first(inc_100k))/first(inc_100k))/n_years),
pct_change_lower_annual = (((last(.lower) - first(.lower))/first(.lower))/n_years),
pct_change_upper_annual = (((last(.upper) - first(.upper))/first(.upper))/n_years),
.by = c(acf_period, {{grouping_var}})) %>%
distinct() %>%
mutate(change = "Slope change")
lst(immediate_change, post_change, slope_change)
}
Function for calculating difference from counterfactual
calcuate_counterfactual <- function(model_data, model, population_denominator, grouping_var=NULL){
#effect vs. counterfactual
counterfact <-
add_epred_draws(object = {{model}},
newdata = {{model_data}} %>%
select(year, year2, y_num, {{population_denominator}}, {{grouping_var}}) %>%
mutate(acf_period = "a. pre-acf")) %>%
group_by(.draw, year, {{grouping_var}}, acf_period) %>%
mutate(.epred_inc_counterf = .epred/{{population_denominator}}*100000, .epred_counterf=.epred) %>%
filter(year>1957) %>%
ungroup() %>%
select(year, {{population_denominator}}, .draw, .epred_counterf, .epred_inc_counterf)
#Calcuate incidence per draw, then summarise.
post_change <-
add_epred_draws(object = {{model}},
newdata = {{model_data}} %>%
select(year, year2, y_num, {{population_denominator}}, {{grouping_var}}, acf_period)) %>%
group_by(.draw, year, {{grouping_var}}, acf_period) %>%
mutate(.epred_inc = .epred/{{population_denominator}}*100000) %>%
filter(year>1957) %>%
ungroup() %>%
select(year, {{population_denominator}}, {{grouping_var}}, .draw, .epred, .epred_inc)
#for the overall period
counterfact_overall <-
add_epred_draws(object = {{model}},
newdata = {{model_data}} %>%
select(year, year2, y_num, {{population_denominator}}, {{grouping_var}}) %>%
mutate(acf_period = "a. pre-acf")) %>%
group_by(.draw, {{grouping_var}}) %>%
filter(year>1957) %>%
ungroup() %>%
select({{population_denominator}}, .draw, .epred) %>%
group_by(.draw) %>%
summarise(.epred_counterf = sum(.epred)) %>%
mutate(year = "Overall (1958-1963)")
#Calcuate incidence per draw, then summarise.
post_change_overall <-
add_epred_draws(object = {{model}},
newdata = {{model_data}} %>%
select(year, year2, y_num, {{population_denominator}}, {{grouping_var}}, acf_period)) %>%
group_by(.draw, {{grouping_var}}) %>%
filter(year>1957) %>%
ungroup() %>%
select({{population_denominator}}, {{grouping_var}}, .draw, .epred) %>%
group_by(.draw, {{grouping_var}}) %>%
summarise(.epred = sum(.epred))
counter_post <-
left_join(counterfact, post_change) %>%
mutate(cases_averted = .epred_counterf-.epred,
pct_change = (.epred - .epred_counterf)/.epred_counterf,
diff_inc100k = .epred_inc - .epred_inc_counterf,
rr_inc100k = .epred_inc/.epred_inc_counterf) %>%
group_by(year, {{grouping_var}}) %>%
mean_qi(cases_averted, pct_change, diff_inc100k, rr_inc100k) %>%
ungroup()
counter_post_overall <-
left_join(counterfact_overall, post_change_overall) %>%
mutate(cases_averted = .epred_counterf-.epred,
pct_change = (.epred - .epred_counterf)/.epred_counterf) %>%
group_by({{grouping_var}}) %>%
mean_qi(cases_averted, pct_change) %>%
ungroup() %>%
mutate(year = "Overall (1958-1963)")
lst(counter_post, counter_post_overall)
}
Function for tidying up counterfactuals (mostly for making nice tables)
tidy_counterfactuals <- function(data){
data %>%
mutate(across(c(cases_averted:cases_averted.upper, diff_inc100k:diff_inc100k.upper), number_format(accuracy = 0.1, big.mark = ","))) %>%
mutate(across(c(rr_inc100k:rr_inc100k.upper), number_format(accuracy = 0.01))) %>%
mutate(across(c(pct_change:pct_change.upper), percent, accuracy=0.1)) %>%
mutate(year = as.character(year),
cases_averted = glue::glue("{cases_averted} ({cases_averted.lower} to {cases_averted.upper})"),
pct_change = glue::glue("{pct_change} ({pct_change.lower} to {pct_change.upper})"),
diff_inc = glue::glue("{diff_inc100k} ({diff_inc100k.lower} to {diff_inc100k.upper})"),
rr_inc = glue::glue("{rr_inc100k} ({rr_inc100k.lower} to {rr_inc100k.upper})"))
}
tidy_counterfactuals_overall <- function(data){
data %>%
mutate(across(c(cases_averted:cases_averted.upper), number_format(accuracy = 0.1, big.mark = ","))) %>%
mutate(across(c(pct_change:pct_change.upper), percent, accuracy=0.1)) %>%
mutate(year = as.character(year),
cases_averted = glue::glue("{cases_averted} ({cases_averted.lower} to {cases_averted.upper})"),
pct_change = glue::glue("{pct_change} ({pct_change.lower} to {pct_change.upper})"))
}
tidy_counterfactuals <- function(data){
data %>%
mutate(across(c(cases_averted:cases_averted.upper, diff_inc100k:diff_inc100k.upper), number_format(accuracy = 0.1, big.mark = ","))) %>%
mutate(across(c(rr_inc100k:rr_inc100k.upper), number_format(accuracy = 0.01))) %>%
mutate(across(c(pct_change:pct_change.upper), percent, accuracy=0.1)) %>%
mutate(year = as.character(year),
cases_averted = glue::glue("{cases_averted} ({cases_averted.lower} to {cases_averted.upper})"),
pct_change = glue::glue("{pct_change} ({pct_change.lower} to {pct_change.upper})"),
diff_inc = glue::glue("{diff_inc100k} ({diff_inc100k.lower} to {diff_inc100k.upper})"),
rr_inc = glue::glue("{rr_inc100k} ({rr_inc100k.lower} to {rr_inc100k.upper})"))
}
Import datasets for analysis
Make a map of Glasgow wards
glasgow_wards_1951 <- st_read(here("mapping/glasgow_wards_1951.geojson"))
Reading layer `glasgow_wards_1951' from data source `/Users/petermacpherson/Dropbox/Projects/Historical TB ACF 2023-11-28/Work/analysis/glasgow-cxr/mapping/glasgow_wards_1951.geojson' using driver `GeoJSON'
Simple feature collection with 37 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -4.393502 ymin: 55.77464 xmax: -4.070411 ymax: 55.92814
Geodetic CRS: WGS 84
#read in Scotland boundary
scotland <- st_read(here("mapping/Scotland_boundary/Scotland boundary.shp"))
Reading layer `Scotland boundary' from data source
`/Users/petermacpherson/Dropbox/Projects/Historical TB ACF 2023-11-28/Work/analysis/glasgow-cxr/mapping/Scotland_boundary/Scotland boundary.shp' using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 1 field
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 5513 ymin: 530249 xmax: 470332 ymax: 1220302
Projected CRS: OSGB36 / British National Grid
#make a bounding box for Glasgow
bbox <- st_bbox(glasgow_wards_1951) |> st_as_sfc()
#plot scotlan with a bounding box around the City of Glasgow
scotland_with_bbox <- ggplot() +
geom_sf(data = scotland, fill="antiquewhite") +
geom_sf(data = bbox, colour = "#C60C30", fill="antiquewhite") +
theme_void() +
theme(panel.border = element_rect(colour = "grey78", fill=NA, linewidth = 0.5),
panel.background = element_rect(fill = "#EAF7FA", size = 0.3))
Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
Please use the `linewidth` argument instead.
#plot the wards
#note we tidy up some names to fit on map
glasgow_ward_map <- glasgow_wards_1951 %>%
mutate(ward = case_when(ward=="Shettleston and Tollcross" ~ "Shettleston and\nTollcross",
ward=="Partick (West)" ~ "Partick\n(West)",
ward=="Partick (East)" ~ "Partick\n(East)",
ward=="North Kelvin" ~ "North\nKelvin",
ward=="Kinning Park" ~ "Kinning\nPark",
TRUE ~ ward)) %>%
ggplot() +
geom_sf(aes(fill=division)) +
geom_sf_label(aes(label = ward), size=3, fill=NA, label.size = NA, colour="black", family = "Segoe UI") +
#scale_colour_identity() +
scale_fill_brewer(palette = "Set3", name="City of Glasgow Division") +
theme_grey(base_family = "Segoe UI") +
labs(x="",
y="",
fill="Division") +
theme(legend.position = "top",
panel.border = element_rect(colour = "grey78", fill=NA, linewidth = 0.5),
panel.background = element_rect(fill = "antiquewhite", size = 0.3),
panel.grid.major = element_line(color = "grey78")) +
guides(fill=guide_legend(title.position = "top", title.hjust = 0.5, title.theme = element_text(face="bold"))) +
scalebar(glasgow_wards_1951, dist = 2, dist_unit = "km",
transform = TRUE, model = "WGS84", location="bottomleft")
#add the map of scotland as an inset
glasgow_ward_map + inset_element(scotland_with_bbox, 0.75, 0, 1, 0.4)
ggsave(here("figures/s1.png"), height=10, width = 12)
NA
NA
Load in the datasets for denonomiators, and check for consistency.
overall_pops <- read_xlsx(path = "2023-11-28_glasgow-acf.xlsx", sheet = "overall_population")
overall_pops %>%
mutate(across(where(is.numeric) & !(year), ~comma(.))) %>%
datatable()
#shift year to midpoint
overall_pops <- overall_pops %>%
mutate(year2 = year+0.5)
Note, we have three population estimates:
(Population in shipping is estimated from the 1951 census, so is the same for most years)
First, plot the total population
overall_pops %>%
ggplot() +
geom_area(aes(y=total_population, x=year2), alpha=0.5, colour = "mediumseagreen", fill="mediumseagreen") +
geom_point(aes(y=total_population, x=year2), colour = "mediumseagreen") +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
scale_x_continuous(labels = year_labels,
breaks = year_labels) +
labs(
title = "Glasgow Corporation: total population",
subtitle = "1950 to 1963",
x = "Year",
y = "Population",
caption = "Mid-year estimates\nMass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)"
) +
theme_ggdist()
NA
NA
Now the population excluding institutionalised and shipping population
overall_pops %>%
ggplot() +
geom_area(aes(y=population_without_inst_ship, x=year2), alpha=0.5, colour = "purple", fill="purple") +
geom_point(aes(y=population_without_inst_ship, x=year2), colour = "purple") +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
scale_x_continuous(labels = year_labels,
breaks = year_labels) +
labs(
title = "Glasgow Corporation: population excluding institutionalised and shipping",
subtitle = "1950 to 1963",
x = "Year",
y = "Population",
caption = "Mid-year estimates\nMass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)"
) +
theme_ggdist()
NA
NA
There are 5 Divisions containing 37 Wards in the Glasgow Corporation, with consistent boundaries over time.
#look-up table for divisions and wards
ward_lookup <- read_xlsx(path = "2023-11-28_glasgow-acf.xlsx", sheet = "divisions_wards")
ward_pops <- read_xlsx(path = "2023-11-28_glasgow-acf.xlsx", sheet = "ward_population")
ward_pops %>%
mutate(across(where(is.numeric) & !(year), ~comma(.))) %>%
datatable()
#shift year to midpoint
ward_pops <- ward_pops %>%
mutate(year2 = year+0.5)
#Get the Division population
division_pops <- ward_pops %>%
group_by(division, year) %>%
summarise(population_without_inst_ship = sum(population_without_inst_ship, na.rm = TRUE),
institutions = sum(institutions, na.rm = TRUE),
shipping = sum(shipping, na.rm = TRUE),
total_population = sum(total_population, na.rm = TRUE))
`summarise()` has grouped output by 'division'. You can override using the `.groups` argument.
division_pops %>%
mutate(across(where(is.numeric) & !(year), ~comma(.))) %>%
datatable()
NA
Plot the overall population by Division and Ward
division_pops %>%
mutate(year2 = year+0.5) %>%
ggplot() +
geom_area(aes(y=total_population, x=year2, colour=division, fill=division), alpha=0.8) +
geom_point(aes(y=total_population, x=year2, colour=division)) +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
facet_wrap(division~.) +
scale_x_continuous(labels = year_labels,
breaks = year_labels,
guide = guide_axis(angle = 90)) +
scale_fill_brewer(palette = "Set3", name = "") +
scale_colour_brewer(palette = "Set3", name = "") +
labs(
title = "Glasgow Corporation: total population by Division",
subtitle = "1950 to 1963",
x = "Year",
y = "Population",
caption = "Mid-year estimates\nMass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)"
) +
theme_ggdist() +
theme(legend.position = "bottom")
NA
NA
ward_pops %>%
ggplot() +
geom_area(aes(y=total_population, x=year2, colour=division, fill=division), alpha=0.5) +
geom_point(aes(y=total_population, x=year2, colour=division)) +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
facet_wrap(ward~., ncol=6) +
scale_x_continuous(labels = year_labels,
breaks = year_labels,
guide = guide_axis(angle = 90)) +
scale_fill_brewer(palette = "Set3", name="Division") +
scale_colour_brewer(palette = "Set3", name = "Division") +
labs(
title = "Glasgow City: total population by Ward",
subtitle = "1950 to 1963",
x = "Year",
y = "Population",
caption = "Mid-year estimates\nMass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)"
) +
theme_ggdist() +
theme(legend.position = "bottom")
ggsave(here("figures/s2.png"), height=10, width=12)
Approximately, how many person-years of follow-up do we have?
overall_pops %>%
ungroup() %>%
summarise(across(year, length, .names = "years"),
across(c(population_without_inst_ship, total_population), sum)) %>%
mutate(across(where(is.double), comma)) %>%
datatable()
NA
NA
Change in population by ward
ward_pops %>%
group_by(ward) %>%
summarise(pct_change_pop = (last(population_without_inst_ship) - first(population_without_inst_ship))/first(population_without_inst_ship)) %>%
mutate(pct_change_pop = percent(pct_change_pop)) %>%
arrange(pct_change_pop) %>%
datatable()
NA
NA
NA
age_sex <- read_xlsx(path = "2023-11-28_glasgow-acf.xlsx", sheet = "age_sex_population") %>%
pivot_longer(cols = c(male, female),
names_to = "sex")
#collapse down to smaller age groups to be manageable
age_sex <- age_sex %>%
ungroup() %>%
mutate(age = case_when(age == "0 to 4" ~ "00 to 04",
age == "5 to 9" ~ "05 to 14",
age == "10 to 14" ~ "05 to 14",
age == "15 to 19" ~ "15 to 24",
age == "20 to 24" ~ "15 to 24",
age == "25 to 29" ~ "25 to 34",
age == "30 to 34" ~ "25 to 34",
age == "35 to 39" ~ "35 to 44",
age == "40 to 44" ~ "35 to 44",
age == "45 to 49" ~ "45 to 59",
age == "50 to 54" ~ "45 to 59",
age == "55 to 59" ~ "45 to 59",
TRUE ~ "60 & up")) %>%
group_by(year, age, sex) %>%
mutate(value = sum(value)) %>%
ungroup()
m_age_sex <- lm(value ~ splines::ns(year, knots = 3)*age*sex, data = age_sex)
summary(m_age_sex)
Warning: essentially perfect fit: summary may be unreliable
Call:
lm(formula = value ~ splines::ns(year, knots = 3) * age * sex,
data = age_sex)
Residuals:
Min 1Q Median 3Q Max
-1.185e-10 0.000e+00 0.000e+00 0.000e+00 1.185e-10
Coefficients: (14 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.222e+04 2.040e-10 2.559e+14 <2e-16 ***
splines::ns(year, knots = 3)1 -8.043e+03 4.071e-10 -1.976e+13 <2e-16 ***
splines::ns(year, knots = 3)2 NA NA NA NA
age05 to 14 3.669e+04 2.499e-10 1.468e+14 <2e-16 ***
age15 to 24 -3.893e+03 2.499e-10 -1.558e+13 <2e-16 ***
age25 to 34 -3.996e+04 2.499e-10 -1.599e+14 <2e-16 ***
age35 to 44 -4.230e+04 2.499e-10 -1.693e+14 <2e-16 ***
age45 to 59 5.459e+04 2.356e-10 2.317e+14 <2e-16 ***
age60 & up 7.533e+04 2.204e-10 3.418e+14 <2e-16 ***
sexmale 3.374e+03 2.886e-10 1.169e+13 <2e-16 ***
splines::ns(year, knots = 3)1:age05 to 14 -1.863e+03 4.985e-10 -3.737e+12 <2e-16 ***
splines::ns(year, knots = 3)2:age05 to 14 NA NA NA NA
splines::ns(year, knots = 3)1:age15 to 24 7.533e+04 4.985e-10 1.511e+14 <2e-16 ***
splines::ns(year, knots = 3)2:age15 to 24 NA NA NA NA
splines::ns(year, knots = 3)1:age25 to 34 1.325e+05 4.985e-10 2.658e+14 <2e-16 ***
splines::ns(year, knots = 3)2:age25 to 34 NA NA NA NA
splines::ns(year, knots = 3)1:age35 to 44 1.380e+05 4.985e-10 2.769e+14 <2e-16 ***
splines::ns(year, knots = 3)2:age35 to 44 NA NA NA NA
splines::ns(year, knots = 3)1:age45 to 59 3.474e+03 4.700e-10 7.390e+12 <2e-16 ***
splines::ns(year, knots = 3)2:age45 to 59 NA NA NA NA
splines::ns(year, knots = 3)1:age60 & up -8.453e+04 4.397e-10 -1.923e+14 <2e-16 ***
splines::ns(year, knots = 3)2:age60 & up NA NA NA NA
splines::ns(year, knots = 3)1:sexmale -1.994e+03 5.757e-10 -3.464e+12 <2e-16 ***
splines::ns(year, knots = 3)2:sexmale NA NA NA NA
age05 to 14:sexmale 1.053e+04 3.534e-10 2.980e+13 <2e-16 ***
age15 to 24:sexmale 2.352e+04 3.534e-10 6.656e+13 <2e-16 ***
age25 to 34:sexmale 1.355e+04 3.534e-10 3.833e+13 <2e-16 ***
age35 to 44:sexmale -1.727e+03 3.534e-10 -4.888e+12 <2e-16 ***
age45 to 59:sexmale 2.774e+03 3.332e-10 8.324e+12 <2e-16 ***
age60 & up:sexmale -7.761e+04 3.117e-10 -2.490e+14 <2e-16 ***
splines::ns(year, knots = 3)1:age05 to 14:sexmale -2.049e+04 7.051e-10 -2.906e+13 <2e-16 ***
splines::ns(year, knots = 3)2:age05 to 14:sexmale NA NA NA NA
splines::ns(year, knots = 3)1:age15 to 24:sexmale -6.780e+04 7.051e-10 -9.616e+13 <2e-16 ***
splines::ns(year, knots = 3)2:age15 to 24:sexmale NA NA NA NA
splines::ns(year, knots = 3)1:age25 to 34:sexmale -3.804e+04 7.051e-10 -5.396e+13 <2e-16 ***
splines::ns(year, knots = 3)2:age25 to 34:sexmale NA NA NA NA
splines::ns(year, knots = 3)1:age35 to 44:sexmale -1.171e+04 7.051e-10 -1.661e+13 <2e-16 ***
splines::ns(year, knots = 3)2:age35 to 44:sexmale NA NA NA NA
splines::ns(year, knots = 3)1:age45 to 59:sexmale -3.473e+04 6.647e-10 -5.224e+13 <2e-16 ***
splines::ns(year, knots = 3)2:age45 to 59:sexmale NA NA NA NA
splines::ns(year, knots = 3)1:age60 & up:sexmale 1.056e+05 6.218e-10 1.698e+14 <2e-16 ***
splines::ns(year, knots = 3)2:age60 & up:sexmale NA NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.074e-11 on 44 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 6.006e+29 on 27 and 44 DF, p-value: < 2.2e-16
age_levels <- age_sex %>% select(age) %>% distinct() %>% pull()
age_sex_nd <-
crossing(
age=age_levels,
sex=c("male", "female"),
year = 1950:1963
)
pred_pops <- age_sex_nd %>% modelr::add_predictions(m_age_sex)
Warning: prediction from a rank-deficient fit may be misleading
pred_pops %>%
ggplot(aes(x=year, y=pred, colour=age)) +
geom_line() +
geom_point() +
facet_grid(sex~.) +
scale_y_continuous(labels = comma, limits = c(0, 125000))
#How well do they match up with our overall populations?
pred_pops %>%
group_by(year) %>%
summarise(sum_pred_pop = sum(pred)) %>%
right_join(overall_pops) %>%
select(year, sum_pred_pop, population_without_inst_ship, total_population) %>%
pivot_longer(cols = c(sum_pred_pop, population_without_inst_ship, total_population)) %>%
ggplot(aes(x=year, y=value, colour=name)) +
geom_point() +
scale_y_continuous(labels = comma, limits = c(800000, 1250000))
Joining with `by = join_by(year)`
pred_pops %>%
group_by(year, sex) %>%
summarise(sum = sum(pred)) %>%
group_by(year) %>%
mutate(sex_ratio = first(sum)/last(sum))
`summarise()` has grouped output by 'year'. You can override using the `.groups` argument.
Population pyramids
label_abs <- function(x) {
comma(abs(x))
}
pred_pops %>%
ungroup() %>%
group_by(year) %>%
mutate(year_pop = sum(pred),
age_sex_pct = percent(pred/year_pop, accuracy=0.1)) %>%
mutate(sex = case_when(sex=="male" ~ "Male",
sex=="female" ~ "Female")) %>%
ggplot(
aes(x = age, fill = sex,
y = ifelse(test = sex == "Female",yes = -pred, no = pred))) +
geom_bar(stat = "identity") +
geom_text(aes(label = age_sex_pct),
position= position_stack(vjust=0.5), colour="white", size=2.5) +
facet_wrap(year~., ncol=7) +
coord_flip() +
scale_y_continuous(labels = label_abs) +
scale_fill_manual(values = c("mediumseagreen", "purple"), name="") +
theme_ggdist() +
theme(axis.text.x = element_text(angle=90, hjust = 1, vjust=0.5),
legend.position = "bottom",
panel.border = element_rect(colour = "grey78", fill=NA)) +
labs(x="", y="")
ggsave(here("figures/s3.png"), width=10)
Saving 10 x 4.51 in image
Not perfect, but resonably good. But ahhhhh… the age groups don’t align with the case notification age groups! Come back to think about this later.
Import the tuberculosis cases dataset
Overall, by year.
cases_by_year <- read_xlsx("2023-11-28_glasgow-acf.xlsx", sheet = "by_year")
cases_by_year%>%
mutate(across(where(is.numeric) & !(year), ~comma(.))) %>%
datatable()
#shift year to midpoint
cases_by_year <- cases_by_year %>%
mutate(year2 = year+0.5)
Plot the overall number of case notified per year, by pulmonary and extra pulmonary classification.
cases_by_year %>%
select(-total_notifications, -year) %>%
pivot_longer(cols = c(pulmonary_notifications, `non-pulmonary_notifications`)) %>%
mutate(name = case_when(name == "pulmonary_notifications" ~ "Pulmonary TB",
name == "non-pulmonary_notifications" ~ "Extra-pulmonary TB")) %>%
ggplot() +
geom_area(aes(y=value, x=year2, group = name, fill=name), alpha=0.5) +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
scale_x_continuous(labels = year_labels,
breaks = year_labels) +
scale_fill_brewer(palette = "Set1", name="") +
labs(
title = "Glasgow Corporation: Tuberculosis notifications",
subtitle = "1950 to 1963, by TB disease classification",
x = "Year",
y = "Number of cases",
caption = "Mass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)"
) +
theme_ggdist() +
theme(legend.position = "bottom")
NA
NA
Read in the datasets and merge together.
#list all the sheets
all_sheets <- excel_sheets("2023-11-28_glasgow-acf.xlsx")
#get the ward sheets
ward_sheets <- enframe(all_sheets) %>%
filter(grepl("by_ward", value)) %>%
pull(value)
cases_by_ward_sex_year <- map_df(ward_sheets, ~read_xlsx(path = "2023-11-28_glasgow-acf.xlsx",
sheet = .))
cases_by_ward_sex_year %>%
mutate(across(where(is.numeric) & !(year), ~comma(.))) %>%
datatable()
NA
Aggregate together to case cases by division
cases_by_division <- cases_by_ward_sex_year %>%
left_join(ward_lookup) %>%
group_by(division, year, tb_type) %>%
summarise(cases = sum(cases, na.rm = TRUE))
Joining with `by = join_by(ward)``summarise()` has grouped output by 'division', 'year'. You can override using the `.groups` argument.
#shift year to midpoint
cases_by_division <- cases_by_division %>%
mutate(year2 = year+0.5) %>%
ungroup()
cases_by_division %>%
select(-year2) %>%
select(year, everything()) %>%
mutate(across(where(is.numeric) & !(year), ~comma(.))) %>%
datatable()
cases_by_division %>%
mutate(tb_type = case_when(tb_type == "Pulmonary" ~ "Pulmonary TB",
tb_type == "Non-Pulmonary" ~ "Extra-pulmonary TB")) %>%
ggplot() +
geom_area(aes(y=cases, x=year2, group = tb_type, fill=tb_type), alpha=0.5) +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
scale_x_continuous(labels = year_labels,
breaks = year_labels,
guide = guide_axis(angle = 90)) +
facet_wrap(division~., scales = "free_y") +
scale_fill_brewer(palette = "Set1", name="") +
labs(
title = "Glasgow Corporation: Tuberculosis notifications by Division",
subtitle = "1950 to 1963, by TB disease classification",
x = "Year",
y = "Number of cases",
caption = "Mass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)\nNote: extra-pulmonary TB cases by Division/Ward not reported in 1962-1963"
) +
theme_ggdist() +
theme(legend.position = "bottom")
cases_by_ward <- cases_by_ward_sex_year %>%
group_by(ward, year, tb_type) %>%
summarise(cases = sum(cases, na.rm = TRUE)) %>%
ungroup()
`summarise()` has grouped output by 'ward', 'year'. You can override using the `.groups` argument.
cases_by_ward %>%
mutate(across(where(is.numeric) & !(year), ~comma(.))) %>%
select(year, everything()) %>%
datatable()
#shift year to midpoint
cases_by_ward <- cases_by_ward %>%
mutate(year2 = year+0.5)
cases_by_ward %>%
mutate(tb_type = case_when(tb_type == "Pulmonary" ~ "Pulmonary TB",
tb_type == "Non-Pulmonary" ~ "Extra-pulmonary TB")) %>%
ggplot() +
geom_area(aes(y=cases, x=year2, group = tb_type, fill=tb_type), alpha=0.8) +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
scale_x_continuous(labels = year_labels,
breaks = year_labels,
guide = guide_axis(angle = 90)) +
facet_wrap(ward~., scales = "free_y") +
scale_fill_brewer(palette = "Set1", name="") +
labs(
title = "Glasgow Corporation: Tuberculosis notifications by Ward",
subtitle = "1950 to 1963, by TB disease classification",
x = "Year",
y = "Number of cases",
caption = "Mass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)\nNote: extra-pulmonary TB cases by Division/Ward not reported in 1962-1963"
) +
theme(legend.position = "bottom")
NA
NA
As we don’t have denominators, we will just model the change in counts.
#list all the sheets
all_sheets <- excel_sheets("2023-11-28_glasgow-acf.xlsx")
#get the ward sheets
age_sex_sheets <- enframe(all_sheets) %>%
filter(grepl("by_age_sex", value)) %>%
pull(value)
cases_by_age_sex <- map_df(age_sex_sheets, ~read_xlsx(path = "2023-11-28_glasgow-acf.xlsx",
sheet = .))
cases_by_age_sex %>%
mutate(across(where(is.numeric) & !(year), ~comma(.))) %>%
datatable()
NA
NA
Now calculate incidence per 100,000 population
Merge the notification and population denominator datasets together.
Here we need to include the whole population (with shipping and institutions) as they are included in the notifications.
overall_inc <- overall_pops %>%
left_join(cases_by_year)
Joining with `by = join_by(year, year2)`
overall_inc <- overall_inc %>%
mutate(inc_pulm_100k = pulmonary_notifications/total_population*100000,
inc_ep_100k = `non-pulmonary_notifications`/total_population*100000,
inc_100k = total_notifications/total_population*100000)
overall_inc %>%
select(year, inc_100k, inc_pulm_100k, inc_ep_100k) %>%
mutate_at(.vars = vars(inc_100k, inc_pulm_100k, inc_ep_100k),
.funs = funs(round)) %>%
datatable()
Warning: `funs()` was deprecated in dplyr 0.8.0.
Please use a list of either functions or lambdas:
# Simple named list:
list(mean = mean, median = median)
# Auto named with `tibble::lst()`:
tibble::lst(mean, median)
# Using lambdas
list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
overall_inc %>%
select(year2, inc_pulm_100k, inc_ep_100k) %>%
pivot_longer(cols = c(inc_pulm_100k, `inc_ep_100k`)) %>%
mutate(name = case_when(name == "inc_pulm_100k" ~ "Pulmonary TB",
name == "inc_ep_100k" ~ "Extra-pulmonary TB")) %>%
ggplot() +
geom_area(aes(y=value, x=year2, group = name, fill=name), alpha=0.5) +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
scale_x_continuous(labels = year_labels,
breaks = year_labels) +
scale_fill_brewer(palette = "Set1", name="") +
labs(
title = "Glasgow Corporation: Tuberculosis case notification rate",
subtitle = "1950 to 1963, by TB disease classification",
x = "Year",
y = "Case notification rate (per 100,000)",
caption = "Mass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)"
) +
theme_ggdist() +
theme(legend.position = "bottom")
NA
NA
NA
division_inc <- division_pops %>%
left_join(cases_by_division)
Joining with `by = join_by(division, year)`
division_inc <- division_inc %>%
mutate(inc_100k = cases/total_population*100000)
division_inc %>%
select(year, division, tb_type, inc_100k) %>%
mutate_at(.vars = vars(inc_100k),
.funs = funs(round)) %>%
datatable()
Warning: `funs()` was deprecated in dplyr 0.8.0.
Please use a list of either functions or lambdas:
# Simple named list:
list(mean = mean, median = median)
# Auto named with `tibble::lst()`:
tibble::lst(mean, median)
# Using lambdas
list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
division_inc %>%
mutate(tb_type = case_when(tb_type == "Pulmonary" ~ "Pulmonary TB",
tb_type == "Non-Pulmonary" ~ "Extra-pulmonary TB")) %>%
ggplot() +
geom_area(aes(y=inc_100k, x=year2, group = tb_type, fill=tb_type), alpha=0.5) +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
scale_x_continuous(labels = year_labels,
breaks = year_labels,
guide = guide_axis(angle = 90)) +
facet_wrap(division~.) +
scale_fill_brewer(palette = "Set1", name="") +
labs(
title = "Glasgow Corporation: Tuberculosis case notification rate, by Division",
subtitle = "1950 to 1963, by TB disease classification",
x = "Year",
y = "Case notification rate (per 100,000)",
caption = "Mass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)\nNote: extra-pulmonary TB cases by Division/Ward not reported in 1962-1963"
) +
theme_ggdist() +
theme(legend.position = "bottom")
NA
NA
NA
Here we will filter out the institutions and harbour from the denominators, as we don’t have reliable population denominators for them.
ward_inc <- ward_pops %>%
left_join(cases_by_ward)
Joining with `by = join_by(ward, year, year2)`
ward_inc <- ward_inc %>%
mutate(inc_100k = cases/population_without_inst_ship*100000)
ward_inc %>%
select(year, ward, tb_type, inc_100k) %>%
mutate_at(.vars = vars(inc_100k),
.funs = funs(round)) %>%
datatable()
Warning: `funs()` was deprecated in dplyr 0.8.0.
Please use a list of either functions or lambdas:
# Simple named list:
list(mean = mean, median = median)
# Auto named with `tibble::lst()`:
tibble::lst(mean, median)
# Using lambdas
list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
ward_inc %>%
mutate(tb_type = case_when(tb_type == "Pulmonary" ~ "Pulmonary TB",
tb_type == "Non-Pulmonary" ~ "Extra-pulmonary TB")) %>%
ggplot() +
geom_area(aes(y=inc_100k, x=year2, group = tb_type, fill=tb_type), alpha=0.5) +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
scale_x_continuous(labels = year_labels,
breaks = year_labels,
guide = guide_axis(angle = 90)) +
facet_wrap(ward~.) +
scale_fill_brewer(palette = "Set1", name="") +
labs(
title = "Glasgow Corporation: Tuberculosis case notification rate, by Ward",
subtitle = "1950 to 1963, by TB disease classification",
x = "Year",
y = "Incidence (per 100,000)",
caption = "Mass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)\nNote: extra-pulmonary TB cases by Division/Ward not reported in 1962-1963"
) +
theme(legend.position = "bottom")
NA
NA
NA
NA
On a map
st_as_sf(left_join(ward_inc, glasgow_wards_1951)) %>%
filter(tb_type=="Pulmonary") %>%
ggplot() +
geom_sf(aes(fill=inc_100k)) +
facet_wrap(year~., ncol = 7) +
scale_fill_viridis_c(name="Case notification rate (per 100,000)",
option = "A") +
theme_ggdist() +
theme(legend.position = "top",
legend.key.width = unit(2, "cm"),
panel.border = element_rect(colour = "grey78", fill=NA)) +
guides(fill=guide_colorbar(title.position = "top"))
Joining with `by = join_by(division, ward, ward_number)`
Import the TB mortality data.
First, overall deaths. Note that in the original reports, we have a pulmonary TB death rate per million for all years, and numbers of pulmonary TB deaths for each year apart from 1950.
#get the overall mortality sheets
deaths_sheets <- enframe(all_sheets) %>%
filter(grepl("deaths", value)) %>%
pull(value)
overall_deaths <- map_df(deaths_sheets, ~read_xlsx(path = "2023-11-28_glasgow-acf.xlsx",
sheet = .))
overall_deaths %>%
mutate(across(where(is.numeric) & !(year), ~comma(.))) %>%
datatable()
NA
NA
NA
Plot the raw numbers of pulmonary deaths
overall_deaths %>%
ggplot(aes(x=year, y=pulmonary_deaths)) +
geom_line(colour = "#DE0D92") +
geom_point(colour = "#DE0D92") +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
labs(y="Pulmonary TB deaths per year",
x = "Year",
title = "Numbers of pulmonary TB deaths",
subtitle = "Glasgow, 1950-1963",
caption = "Mass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)\nNote: no data for 1950") +
theme_ggdist() +
theme(panel.border = element_rect(colour = "grey78", fill=NA))
NA
NA
Now the incidence of pulmonary TB death
overall_deaths %>%
ggplot(aes(x=year, y=pulmonary_death_rate_per_100k)) +
geom_line(colour = "#4D6CFA") +
geom_point(colour = "#4D6CFA") +
geom_vline(aes(xintercept=acf_start), linetype=3) +
geom_vline(aes(xintercept=acf_end), linetype=3) +
scale_y_continuous(labels=comma) +
scale_x_continuous(labels = year_labels,
breaks = year_labels) +
labs(y="Annual incidence of death (per 100,000)",
x = "Year",
caption = "Mass miniature X-ray campaign period between dashed lines (11th March-12th April 1957)") +
theme_ggdist() +
theme(panel.border = element_rect(colour = "grey78", fill=NA))
ggsave(here("figures/s7.png"), width=10)
Saving 10 x 4.51 in image
Make Table 1 here, and save for publication.
overall_pops %>%
select(year, total_population) %>%
left_join(overall_inc %>%
select(year,
pulmonary_notifications, inc_pulm_100k,
`non-pulmonary_notifications`, inc_ep_100k,
total_notifications, inc_100k)) %>%
left_join(overall_deaths %>%
select(year,
pulmonary_deaths, pulmonary_death_rate_per_100k)) %>%
mutate(across(where(is.numeric) & !(year), ~round(., digits=1))) %>%
mutate(across(where(is.numeric) & !(year), ~comma(.)))
Joining with `by = join_by(year)`Joining with `by = join_by(year)`
First model will investigate the impact of mass miniature X-ray campaign on pulmonary TB case notification rate using an interrupted time series analysis.
Set up the data
mdata1 <- overall_inc %>%
mutate(acf_period = case_when(year %in% c(1950:1956) ~ "a. pre-acf",
year %in% c(1957) ~ "b. acf",
year %in% c(1958:1963) ~ "c. post-acf")) %>%
mutate(y_num = 1:nrow(.)) %>%
rename(extrapulmonary_notifications = `non-pulmonary_notifications`)
Work on the priors a bit
Basic prior
basic_prior <- c(prior(normal(0, 1), class = Intercept),
prior(normal(0, 0.25), class = b))
Look at the mean and variance of counts (counts of pulmonary notifications are what we are predicting)
#Mean of counts per year
mean(mdata1$pulmonary_notifications)
[1] 1858.429
#variance of counts per year
var(mdata1$pulmonary_notifications)
[1] 716579.8
Quite a bit of over-dispersion here, so negative binomial distribution might be a better choice of distributional family than Poisson.
Slightly more informative prior (“weakly informative” really)
ggplot(data = tibble(x = seq(from = 500, to = 5000, by = 10)),
aes(x = x, y = dgamma(x, shape = 0.1, rate = 0.00001))) +
geom_area(color = "transparent",
fill = "#DE0D92") +
scale_x_continuous(NULL) +
coord_cartesian(xlim = c(500, 5000)) +
ggtitle(expression(brms~~gamma(0.1*", "*0.00001)~shape~prior))
NA
NA
Fit a model with only priors
winform_prior <- c(prior(normal(0, 1), class = Intercept),
prior(gamma(0.1, 0.00001), class = shape),
prior(normal(0, 0.25), class = b))
m_pulmonary_prior <- brm(
pulmonary_notifications ~ y_num + acf_period + acf_period:y_num + offset(log(total_population)),
data = mdata1,
family = negbinomial(),
seed = 1234,
chains = 4, cores = 4,
prior = winform_prior,
sample_prior = "only",
backend="cmdstanr",
save_pars = save_pars(all = TRUE),
warmup = 1000)
Compiling Stan program...
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
In file included from /var/folders/bh/wr0g5x9j2wq46p9hm3ft_b380000gn/T/RtmpvqbQUj/model-115da3e1ad624.hpp:1:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/src/stan/model/model_header.hpp:4:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math.hpp:19:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/rev.hpp:10:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/rev/fun.hpp:198:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor.hpp:15:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor/integrate_ode_rk45.hpp:6:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:9:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint.hpp:76:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/observer_collection.hpp:23:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function.hpp:30:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function/detail/prologue.hpp:17:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function/function_base.hpp:21:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/type_index.hpp:29:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/type_index/stl_type_index.hpp:47:
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.
/
0/boost/container_hash/hash.hpp:132:33: warning: 'unary_function<const std::error_category *, unsigned long>' is deprecated [-Wdeprecated-declarations]
-
struct hash_base : std::unary_function<T, std::size_t> {};
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:692:18: note: in instantiation of template class 'boost::hash_detail::hash_base<const std::error_category *>' requested here
: public boost::hash_detail::hash_base<T*>
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:420:24: note: in instantiation of template class 'boost::hash<const std::error_category *>' requested here
boost::hash<T> hasher;
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:551:9: note: in instantiation of function template specialization 'boost::hash_combine<const std::error_category *>' requested here
hash_combine(seed, &v.category());
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__functional/unary_function.h:23:29: note: 'unary_function<const std::error_category *, unsigned long>' has been explicitly marked deprecated here
struct _LIBCPP_TEMPLATE_VIS _LIBCPP_DEPRECATED_IN_CXX11 unary_function
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__config:850:41: note: expanded from macro '_LIBCPP_DEPRECATED_IN_CXX11'
# define _LIBCPP_DEPRECATED_IN_CXX11 _LIBCPP_DEPRECATED
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__config:835:49: note: expanded from macro '_LIBCPP_DEPRECATED'
# define _LIBCPP_DEPRECATED __attribute__((deprecated))
^
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
1 warning generated.
/
-
ld: warning: duplicate -rpath '/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
\
|
/
Start sampling
Running MCMC with 4 parallel chains...
Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.0 seconds.
Chain 3 finished in 0.1 seconds.
Chain 4 finished in 0.1 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.6 seconds.
Warning: 86 of 4000 (2.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.
summary(m_pulmonary_prior)
Warning: There were 86 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Family: negbinomial
Links: mu = log; shape = identity
Formula: pulmonary_notifications ~ y_num + acf_period + acf_period:y_num + offset(log(total_population))
Data: mdata1 (Number of observations: 14)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.05 2.48 -5.09 4.62 1.00 2538 2031
y_num 0.00 0.25 -0.50 0.48 1.00 2589 2402
acf_periodb.acf 0.00 0.25 -0.49 0.48 1.00 2908 1869
acf_periodc.postMacf 0.01 0.25 -0.48 0.49 1.00 2668 2293
y_num:acf_periodb.acf -0.00 0.25 -0.49 0.46 1.00 2714 2401
y_num:acf_periodc.postMacf 0.00 0.26 -0.51 0.54 1.00 2137 1486
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 11729.33 36453.14 0.00 104278.62 1.01 937 801
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
conditional_effects(m_pulmonary_prior)
Now fit the model with the weakly informative priors
m_pulmonary_overall <- brm(
pulmonary_notifications ~ y_num + acf_period + acf_period:y_num + offset(log(total_population)),
data = mdata1,
family = negbinomial(),
seed = 1234,
chains = 4, cores = 4,
prior = winform_prior,
save_pars = save_pars(all = TRUE),
backend = "cmdstanr",
warmup = 1000)
Compiling Stan program...
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
In file included from /var/folders/bh/wr0g5x9j2wq46p9hm3ft_b380000gn/T/RtmpvqbQUj/model-115da2adb5f52.hpp:1:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/src/stan/model/model_header.hpp:4:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math.hpp:19:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/rev.hpp:10:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/rev/fun.hpp:198:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor.hpp:15:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor/integrate_ode_rk45.hpp:6:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:9:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint.hpp:76:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/observer_collection.hpp:23:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function.hpp:30:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function/detail/prologue.hpp:17:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function/function_base.hpp:21:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/type_index.hpp:29:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/type_index/stl_type_index.hpp:47:
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.
-
0/boost/container_hash/hash.hpp:132:33: warning: 'unary_function<const std::error_category *, unsigned long>' is deprecated [-Wdeprecated-declarations]
struct hash_base : std::unary_function<T, std::size_t> {};
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:692:18: note: in instantiation of template class 'boost::hash_detail::hash_base<const std::error_category *>' requested here
: public boost::hash_detail::hash_base<T*>
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:420:24: note: in instantiation of template class 'boost::hash<const std::error_category *>' requested here
boost::hash<T> hasher;
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:551:9: note: in instantiation of function template specialization 'boost::hash_combine<const std::error_category *>' requested here
hash_combine(seed, &v.category());
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__functional/unary_function.h:23:29: note: 'unary_function<const std::error_category *, unsigned long>' has been explicitly marked deprecated here
\
struct _LIBCPP_TEMPLATE_VIS _LIBCPP_DEPRECATED_IN_CXX11 unary_function
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__config:850:41: note: expanded from macro '_LIBCPP_DEPRECATED_IN_CXX11'
# define _LIBCPP_DEPRECATED_IN_CXX11 _LIBCPP_DEPRECATED
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__config:835:49: note: expanded from macro '_LIBCPP_DEPRECATED'
# define _LIBCPP_DEPRECATED __attribute__((deprecated))
^
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
1 warning generated.
/
-
ld: warning: duplicate -rpath '/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
\
|
/
Start sampling
Running MCMC with 4 parallel chains...
Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 1 finished in 0.5 seconds.
Chain 2 finished in 0.4 seconds.
Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 finished in 0.4 seconds.
Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 finished in 0.5 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.5 seconds.
Total execution time: 0.8 seconds.
summary(m_pulmonary_overall)
Family: negbinomial
Links: mu = log; shape = identity
Formula: pulmonary_notifications ~ y_num + acf_period + acf_period:y_num + offset(log(total_population))
Data: mdata1 (Number of observations: 14)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -6.10 0.03 -6.17 -6.04 1.00 2999 2234
y_num -0.02 0.01 -0.03 -0.01 1.00 2904 2266
acf_periodb.acf 0.02 0.25 -0.45 0.53 1.00 1764 2081
acf_periodc.postMacf 0.05 0.11 -0.19 0.26 1.00 1935 1892
y_num:acf_periodb.acf 0.08 0.03 0.01 0.14 1.00 1755 2160
y_num:acf_periodc.postMacf -0.05 0.01 -0.08 -0.03 1.00 1821 1648
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 4409.62 15136.27 322.80 26935.74 1.00 1318 1213
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(m_pulmonary_overall)
pp_check(m_pulmonary_overall, type='ecdf_overlay')
Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
Summarise the posterior
f1b <- plot_counterfactual(model_data = mdata1, model = m_pulmonary_overall, population_denominator = total_population, outcome = inc_pulm_100k, grouping_var=NULL)
f1b
Make this into a figure
f1a <- st_as_sf(left_join(ward_inc, glasgow_wards_1951)) %>%
filter(tb_type=="Pulmonary") %>%
ggplot() +
geom_sf(aes(fill=inc_100k)) +
facet_wrap(year~., ncol = 7) +
scale_fill_viridis_c(name="Case notification rate (per 100,000)",
option = "A") +
theme_ggdist() +
theme(panel.border = element_rect(colour = "grey78", fill=NA),
legend.position = "top",
#legend.key.width = unit(2, "cm"),
legend.title.align = 0.5) +
guides(fill=guide_colorbar(title.position = "top"))
Joining with `by = join_by(division, ward, ward_number)`
(f1a / f1b) + plot_annotation(tag_levels = "A")
ggsave(here("figures/f1.png"))
Saving 7.29 x 4.51 in image
Summary of change in notifications
summarise_change(model_data=mdata1, model=m_pulmonary_overall, population_denominator=total_population, grouping_var=NULL) %>%
map(datatable)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.
$immediate_change
$post_change
$slope_change
NA
(Alternative way - keep in for now)
overall_immediate_draws <- mdata1 %>%
select(year, year2, y_num, acf_period, total_population, pulmonary_notifications) %>%
filter(year %in% c(1956,1957)) %>%
add_epred_draws(m_pulmonary_overall) %>%
mutate(inc_100k = .epred/total_population*100000) %>%
group_by(.draw) %>%
summarise(pct_change_immediate = (last(inc_100k) - first(inc_100k))/first(inc_100k)) %>%
ungroup()
overall_post_draws <- mdata1 %>%
select(year, year2, y_num, acf_period, total_population, pulmonary_notifications) %>%
filter(year %in% c(1956,1958)) %>%
add_epred_draws(m_pulmonary_overall) %>%
mutate(inc_100k = .epred/total_population*100000) %>%
group_by(.draw) %>%
summarise(pct_change_post = (last(inc_100k) - first(inc_100k))/first(inc_100k)) %>%
ungroup()
overall_slope_draws <- mdata1 %>%
select(year, year2, y_num, acf_period, total_population, pulmonary_notifications) %>%
filter(year %in% c(1950, 1956, 1958, 1963)) %>%
add_epred_draws(m_pulmonary_overall) %>%
mutate(inc_100k = .epred/total_population*100000) %>%
ungroup() %>%
mutate(n_years = length(year), .by=acf_period) %>%
group_by(acf_period, .draw) %>%
summarise(pct_change_slope = ((last(inc_100k) - first(inc_100k))/first(inc_100k))/n_years) %>%
distinct() %>%
pivot_wider(names_from = c(acf_period),
values_from = pct_change_slope) %>%
mutate(ratio_annual_slope = `c. post-acf` / `a. pre-acf`)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.`summarise()` has grouped output by 'acf_period', '.draw'. You can override using the `.groups` argument.
Correlation between immediate effect and post effect of ACF
left_join(overall_immediate_draws, overall_post_draws) %>%
ggplot(aes(x=pct_change_immediate, y=pct_change_post)) +
geom_hdr(
aes(fill = after_stat(probs)),
alpha = 1) +
#geom_hdr_points(aes(colour = after_stat(probs)), size=0.5) +
geom_smooth(method = "lm", se=FALSE, colour="black", linewidth=0.5) +
stat_regline_equation(label.x = 0.25, label.y = -0.25, size=4) +
scale_x_continuous(labels = percent,
breaks = pretty_breaks(n = 10)) +
scale_y_continuous(labels = percent,
breaks = pretty_breaks(n = 5)) +
scale_fill_viridis_d(option="E", name="") +
labs(title="Correlation between immediate ACF impact and post-ACF case notification rate",
y="Post intervention impact: ercentage change in CNR (1958 vs. 1956)",
x="Immediate impact: percentage change in CNR (1957 vs. 1956)",
caption="Boundaries are posterior desnity intervals from 4000 draws") +
theme_ggdist() +
theme(legend.position = "bottom",
panel.border = element_rect(colour = "grey78", fill=NA))
Joining with `by = join_by(.draw)`
Correlation between immediate effect and change in slope
left_join(overall_immediate_draws, overall_slope_draws) %>%
ggplot(aes(x=pct_change_immediate, y=ratio_annual_slope)) +
geom_hdr(
aes(fill = after_stat(probs)),
alpha = 1) +
#geom_hdr_points(aes(colour = after_stat(probs)), size=0.5) +
geom_smooth(method = "lm", se=FALSE, colour="black", linewidth=0.5) +
#stat_regline_equation(label.x = 0.25, label.y = 0.02, size=4) +
scale_x_continuous(labels = percent,
breaks = pretty_breaks(n = 10)) +
scale_y_continuous(breaks = pretty_breaks(n = 5),
limits = c(0, 10)) +
scale_fill_viridis_d(option="E", name="") +
labs(title="Correlation between immediate ACF impact and post-ACF case notification rate",
y="Post intervention impact: Percentage change in CNR (1958 vs. 1956)",
x="Immediate impact: percentage change in CNR (1957 vs. 1956)",
caption="Points are draws from posteior distribution") +
theme_ggdist() +
theme(legend.position = "bottom",
panel.border = element_rect(colour = "grey78", fill=NA))
Joining with `by = join_by(.draw)`
overall_pulmonary_counterf <- calcuate_counterfactual(model_data = mdata1, model=m_pulmonary_overall, population_denominator = total_population)
Joining with `by = join_by(year, total_population, .draw)`Joining with `by = join_by(.draw)`
overall_pulmonary_counterf$counter_post %>%
mutate(across(c(cases_averted:cases_averted.upper, diff_inc100k:diff_inc100k.upper), number_format(accuracy = 0.1, big.mark = ","))) %>%
mutate(across(c(rr_inc100k:rr_inc100k.upper), number_format(accuracy = 0.01))) %>%
mutate(across(c(pct_change:pct_change.upper), percent, accuracy=0.1)) %>%
datatable()
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `across(c(pct_change:pct_change.upper), percent, accuracy = 0.1)`.
Caused by warning:
! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
Supply arguments directly to `.fns` through an anonymous function instead.
# Previously
across(a:b, mean, na.rm = TRUE)
# Now
across(a:b, \(x) mean(x, na.rm = TRUE))
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
Total pulmonary TB cases averted between 1958 and 1963
overall_pulmonary_counterf$counter_post_overall %>%
mutate(across(c(cases_averted:cases_averted.upper), number_format(accuracy = 0.1, big.mark = ","))) %>%
mutate(across(c(pct_change:pct_change.upper), percent, accuracy=0.1)) %>%
datatable()
NA
NA
m_extrap_overall <- brm(
extrapulmonary_notifications ~ y_num + acf_period + acf_period:y_num + offset(log(total_population)),
data = mdata1,
family = poisson(),
seed = 1234,
chains = 4, cores = 4,
prior = basic_prior,
save_pars = save_pars(all = TRUE),
backend = "cmdstanr",
warmup = 1000)
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
Running MCMC with 4 parallel chains...
Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1 finished in 0.2 seconds.
Chain 2 finished in 0.2 seconds.
Chain 3 finished in 0.2 seconds.
Chain 4 finished in 0.2 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.2 seconds.
Total execution time: 0.5 seconds.
summary(m_extrap_overall)
Family: poisson
Links: mu = log
Formula: extrapulmonary_notifications ~ y_num + acf_period + acf_period:y_num + offset(log(total_population))
Data: mdata1 (Number of observations: 14)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -7.90 0.05 -7.99 -7.81 1.00 2908 2532
y_num -0.09 0.01 -0.11 -0.07 1.00 2623 2183
acf_periodb.acf -0.00 0.24 -0.48 0.48 1.00 2558 2487
acf_periodc.postMacf -0.32 0.17 -0.64 0.01 1.00 2095 2415
y_num:acf_periodb.acf -0.02 0.03 -0.08 0.04 1.00 2413 2577
y_num:acf_periodc.postMacf 0.02 0.02 -0.02 0.05 1.00 1711 2108
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(m_extrap_overall)
pp_check(m_extrap_overall, type='ecdf_overlay')
plot_counterfactual(model_data = mdata1, model=m_extrap_overall, population_denominator = total_population, outcome=inc_ep_100k)
ggsave(here("figures/s6.png"), width=10)
Saving 10 x 4.51 in image
A. Percentage change in mortality, from 1956 to 1957 (i.e. immediate ACF effect)
summarise_change(model_data=mdata1, model = m_extrap_overall, population_denominator = total_population) %>%
map(datatable)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.
$immediate_change
$post_change
$slope_change
NA
overall_ep_counterf <- calcuate_counterfactual(model_data = mdata1, model=m_extrap_overall, population_denominator = total_population)
Joining with `by = join_by(year, total_population, .draw)`Joining with `by = join_by(.draw)`
overall_ep_counterf$counter_post %>%
mutate(across(c(cases_averted:cases_averted.upper, diff_inc100k:diff_inc100k.upper), number_format(accuracy = 0.1, big.mark = ","))) %>%
mutate(across(c(rr_inc100k:rr_inc100k.upper), number_format(accuracy = 0.01))) %>%
mutate(across(c(pct_change:pct_change.upper), percent, accuracy=0.1)) %>%
datatable()
NA
Total extra pulmonary TB cases averted between 1958 and 1963
overall_ep_counterf$counter_post_overall %>%
mutate(across(c(cases_averted:cases_averted.upper), number_format(accuracy = 0.1, big.mark = ","))) %>%
mutate(across(c(pct_change:pct_change.upper), percent, accuracy=0.1)) %>%
datatable()
NA
NA
mdata2 <- ward_inc %>%
filter(tb_type=="Pulmonary") %>%
mutate(acf_period = case_when(year %in% c(1950:1956) ~ "a. pre-acf",
year %in% c(1957) ~ "b. acf",
year %in% c(1958:1963) ~ "c. post-acf")) %>%
group_by(ward) %>%
mutate(y_num = row_number()) %>%
ungroup()
(Note the denominator without institutionalised people and “shipping”!)
#weakly informative priors for multilevel model
basic_prior2 <- c(prior(normal(0, 1), class = Intercept),
prior(normal(0, 0.1), class = b),
prior(cauchy(0,5), class="sd"))
ggplot(data = tibble(x = seq(from = 0, to = 250, by = 10)),
aes(x = x, y = dgamma(x, shape = 0.5, rate = 0.01))) +
geom_area(color = "transparent",
fill = "#DE0D92") +
scale_x_continuous(NULL) +
scale_y_continuous(NULL, breaks = NULL) +
coord_cartesian(xlim = c(0, 250)) +
ggtitle(expression(brms~~gamma(0.5*", "*0.00001)~shape~prior))
winform_prior2 <- c(prior(normal(0, 1), class = Intercept),
prior(gamma(0.5, 0.01), class = shape),
prior(normal(0, 1), class = b),
prior(cauchy(0,5), class="sd"),
prior(lkj(2), class="cor"))
m_pulmonary_ward_prior <- brm(
cases ~ y_num + acf_period + acf_period:y_num + (1 + y_num + acf_period + acf_period:y_num | ward) + offset(log(population_without_inst_ship)),
data = mdata2,
family = negbinomial(),
seed = 1234,
chains = 4, cores = 4,
prior = winform_prior2,
save_pars = save_pars(all = TRUE),
sample_prior = "only",
backend = "cmdstanr",
warmup = 1000)
Compiling Stan program...
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
In file included from /var/folders/bh/wr0g5x9j2wq46p9hm3ft_b380000gn/T/RtmpvqbQUj/model-115da63cc2d2b.hpp:1:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/src/stan/model/model_header.hpp:4:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math.hpp:19:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/rev.hpp:10:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/rev/fun.hpp:198:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor.hpp:15:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor/integrate_ode_rk45.hpp:6:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:9:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint.hpp:76:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/observer_collection.hpp:23:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function.hpp:30:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function/detail/prologue.hpp:17:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function/function_base.hpp:21:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/type_index.hpp:29:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/type_index/stl_type_index.hpp:47:
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.
-
0/boost/container_hash/hash.hpp:132:33: warning: 'unary_function<const std::error_category *, unsigned long>' is deprecated [-Wdeprecated-declarations]
struct hash_base : std::unary_function<T, std::size_t> {};
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:692:18: note: in instantiation of template class 'boost::hash_detail::hash_base<const std::error_category *>' requested here
: public boost::hash_detail::hash_base<T*>
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:420:24: note: in instantiation of template class 'boost::hash<const std::error_category *>' requested here
boost::hash<T> hasher;
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:551:9: note: in instantiation of function template specialization 'boost::hash_combine<const std::error_category *>' requested here
hash_combine(seed, &v.category());
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__functional/unary_function.h:23:29: note: 'unary_function<const std::error_category *, unsigned long>' has been explicitly marked deprecated here
struct _LIBCPP_TEMPLATE_VIS _LIBCPP_DEPRECATED_IN_CXX11 unary_function
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__config:850:41: note: expanded from macro '_LIBCPP_DEPRECATED_IN_CXX11'
# define _LIBCPP_DEPRECATED_IN_CXX11 _LIBCPP_DEPRECATED
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__config:835:49: note: expanded from macro '_LIBCPP_DEPRECATED'
# define _LIBCPP_DEPRECATED __attribute__((deprecated))
^
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
1 warning generated.
|
/
ld: warning: duplicate -rpath '/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
-
\
|
Start sampling
Running MCMC with 4 parallel chains...
Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1 finished in 0.8 seconds.
Chain 2 finished in 0.8 seconds.
Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 finished in 0.9 seconds.
Chain 4 finished in 0.9 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.8 seconds.
Total execution time: 1.1 seconds.
conditional_effects(m_pulmonary_ward_prior)
NA
NA
Now fit the model with data
m_pulmonary_ward <- brm(
cases ~ y_num + acf_period + acf_period:y_num + (1 + y_num + acf_period + acf_period:y_num | ward) + offset(log(population_without_inst_ship)),
data = mdata2,
family = negbinomial(),
seed = 1234,
chains = 4, cores = 4,
prior = winform_prior2,
backend = "cmdstanr")
Compiling Stan program...
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
In file included from /var/folders/bh/wr0g5x9j2wq46p9hm3ft_b380000gn/T/RtmpvqbQUj/model-115da516879bf.hpp:1:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/src/stan/model/model_header.hpp:4:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math.hpp:19:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/rev.hpp:10:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/rev/fun.hpp:198:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor.hpp:15:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor/integrate_ode_rk45.hpp:6:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:9:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint.hpp:76:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/observer_collection.hpp:23:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function.hpp:30:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function/detail/prologue.hpp:17:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function/function_base.hpp:21:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/type_index.hpp:29:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/type_index/stl_type_index.hpp:47:
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.
-
0/boost/container_hash/hash.hpp:132:33: warning: 'unary_function<const std::error_category *, unsigned long>' is deprecated [-Wdeprecated-declarations]
struct hash_base : std::unary_function<T, std::size_t> {};
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:692:18: note: in instantiation of template class 'boost::hash_detail::hash_base<const std::error_category *>' requested here
: public boost::hash_detail::hash_base<T*>
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:420:24: note: in instantiation of template class 'boost::hash<const std::error_category *>' requested here
boost::hash<T> hasher;
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:551:9: note: in instantiation of function template specialization 'boost::hash_combine<const std::error_category *>' requested here
hash_combine(seed, &v.category());
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__functional/unary_function.h:23:29: note: 'unary_function<const std::error_category *, unsigned long>' has been explicitly marked deprecated here
struct _LIBCPP_TEMPLATE_VIS _LIBCPP_DEPRECATED_IN_CXX11 unary_function
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__config:850:41: note: expanded from macro '_LIBCPP_DEPRECATED_IN_CXX11'
\
# define _LIBCPP_DEPRECATED_IN_CXX11 _LIBCPP_DEPRECATED
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__config:835:49: note: expanded from macro '_LIBCPP_DEPRECATED'
# define _LIBCPP_DEPRECATED __attribute__((deprecated))
^
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
1 warning generated.
/
-
ld: warning: duplicate -rpath '/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
\
|
/
Start sampling
Running MCMC with 4 parallel chains...
Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2 finished in 68.7 seconds.
Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 finished in 69.8 seconds.
Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 finished in 70.2 seconds.
Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1 finished in 71.8 seconds.
All 4 chains finished successfully.
Mean chain execution time: 70.1 seconds.
Total execution time: 71.9 seconds.
summary(m_pulmonary_ward)
Family: negbinomial
Links: mu = log; shape = identity
Formula: cases ~ y_num + acf_period + acf_period:y_num + (1 + y_num + acf_period + acf_period:y_num | ward) + offset(log(population_without_inst_ship))
Data: mdata2 (Number of observations: 518)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Group-Level Effects:
~ward (Number of levels: 37)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.26 0.04 0.20 0.34 1.00 1698 2316
sd(y_num) 0.02 0.01 0.01 0.03 1.01 1217 978
sd(acf_periodb.acf) 0.08 0.05 0.00 0.18 1.00 1573 2064
sd(acf_periodc.postMacf) 0.13 0.08 0.01 0.30 1.00 939 1471
sd(y_num:acf_periodb.acf) 0.01 0.01 0.00 0.02 1.00 1612 2176
sd(y_num:acf_periodc.postMacf) 0.01 0.01 0.00 0.03 1.01 576 1121
cor(Intercept,y_num) -0.51 0.20 -0.81 -0.04 1.00 2566 2373
cor(Intercept,acf_periodb.acf) -0.28 0.33 -0.80 0.43 1.00 2987 2293
cor(y_num,acf_periodb.acf) -0.06 0.32 -0.64 0.58 1.00 5536 3410
cor(Intercept,acf_periodc.postMacf) -0.18 0.28 -0.67 0.41 1.00 4468 3115
cor(y_num,acf_periodc.postMacf) 0.14 0.30 -0.45 0.70 1.00 3192 3324
cor(acf_periodb.acf,acf_periodc.postMacf) 0.10 0.32 -0.54 0.69 1.00 1909 2968
cor(Intercept,y_num:acf_periodb.acf) -0.27 0.33 -0.79 0.45 1.00 3493 2858
cor(y_num,y_num:acf_periodb.acf) -0.07 0.32 -0.66 0.57 1.00 5401 3217
cor(acf_periodb.acf,y_num:acf_periodb.acf) -0.09 0.34 -0.71 0.56 1.00 3930 3316
cor(acf_periodc.postMacf,y_num:acf_periodb.acf) 0.07 0.33 -0.58 0.67 1.00 3674 3348
cor(Intercept,y_num:acf_periodc.postMacf) 0.02 0.31 -0.58 0.61 1.00 4641 2858
cor(y_num,y_num:acf_periodc.postMacf) -0.10 0.33 -0.69 0.57 1.00 1930 2364
cor(acf_periodb.acf,y_num:acf_periodc.postMacf) 0.07 0.33 -0.58 0.65 1.00 2736 3112
cor(acf_periodc.postMacf,y_num:acf_periodc.postMacf) -0.16 0.36 -0.81 0.55 1.00 1594 1728
cor(y_num:acf_periodb.acf,y_num:acf_periodc.postMacf) 0.07 0.33 -0.57 0.68 1.00 1848 2837
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -6.14 0.05 -6.24 -6.04 1.00 1155 1961
y_num -0.02 0.01 -0.03 -0.01 1.00 3100 3310
acf_periodb.acf 0.02 1.02 -1.99 1.98 1.00 4602 3119
acf_periodc.postMacf 0.05 0.11 -0.16 0.25 1.00 4476 2977
y_num:acf_periodb.acf 0.08 0.13 -0.16 0.33 1.00 4578 3148
y_num:acf_periodc.postMacf -0.05 0.01 -0.07 -0.03 1.00 4216 3135
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 98.15 23.53 63.84 158.13 1.00 3107 2818
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(m_pulmonary_ward)
pp_check(m_pulmonary_ward, type='ecdf_overlay')
Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
plot_counterfactual(model_data = mdata2, model=m_pulmonary_ward, outcome = inc_100k, population_denominator = population_without_inst_ship, grouping_var = ward, ward)
ggsave(here("figures/s4.png"), width=10, height=12)
A. percentage increase in CNR, from 1956 to 1957 (i.e. immediate ACF effect)
ward_change <- summarise_change(model_data = mdata2, model = m_pulmonary_ward, population_denominator = population_without_inst_ship, grouping_var=ward)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.
ward_change %>%
map(datatable)
$immediate_change
$post_change
$slope_change
NA
As a supplementary figure
ward_change$immediate_change %>%
arrange(acf_inc100k_rr) %>%
ggplot() +
geom_pointrange(aes(y=acf_inc100k_rr, ymin=acf_inc100k_rr.lower, ymax=acf_inc100k_rr.upper,
x=fct_reorder(ward, acf_inc100k_rr),
colour = acf_inc100k_rr)) +
geom_hline(aes(yintercept=1), linetype=2) +
coord_flip() +
scale_colour_viridis_b(option = "D") +
scale_y_continuous(limits = c(0.8,3.0)) +
labs(x="",
y="Relative posterior predicted case notification rate (per 100,000; 95% UI)\nACF (1957) vs. Before ACF (1956)") +
theme_ggdist() +
theme(legend.position = "none",
panel.border = element_rect(colour = "grey78", fill=NA))
ggsave(here("figures/s5.png"))
Saving 7.29 x 4.51 in image
ward_change$post_change %>%
arrange(acf_inc100k_rr) %>%
ggplot() +
geom_pointrange(aes(y=acf_inc100k_rr, ymin=acf_inc100k_rr.lower, ymax=acf_inc100k_rr.upper,
x=fct_reorder(ward, acf_inc100k_rr),
colour = acf_inc100k_rr)) +
geom_hline(aes(yintercept=1), linetype=2) +
coord_flip() +
scale_colour_viridis_b(option = "D") +
#scale_y_continuous(limits = c(0.8,3.0)) +
labs(x="",
y="Relative posterior predicted case notification rate (per 100,000; 95% UI)\nAfter ACF (1958) vs. Before ACF (1956)") +
theme_ggdist() +
theme(legend.position = "none",
panel.border = element_rect(colour = "grey78", fill=NA))
ggsave(here("figures/s6.png"))
Saving 7.29 x 4.51 in image
percentage change = (final value - initial value) / initial value
ward_change$slope_change %>%
ggplot() +
geom_pointrange(aes(y=pct_change_epred_overall , ymin=pct_change_lower_overall , ymax=pct_change_upper_overall ,
group=acf_period, colour=acf_period,
x = fct_reorder(ward, pct_change_epred_overall ))) +
coord_flip() +
scale_y_continuous(labels =percent) +
scale_colour_manual(values = c("#DE0D92", "#4D6CFA")) +
#scale_y_continuous(limits = c(0.8,3.0)) +
labs(title="Intervention effect of mass miniature chest X-ray in Glasgow",
subtitle="By municipal ward",
x="",
y="Mean annual rate of change in case notification rate (95% CrI)\n Before ACF (1950-1956) vs. after ACF (1958-1963)") +
theme_ggdist() +
theme(panel.border = element_rect(colour = "grey78", fill=NA))
(Alternative figure - keep in for the minute)
Is there any correlation between immediate increase and a) post-intervention (1958) effect, and b) post intervention slope (1958-1963)
ward_cors <- ward_impact_out %>%
select(ward, immediate_effect = acf_inc100k_rr,
immediate_effect_lower = acf_inc100k_rr.lower,
immediate_effect_upper = acf_inc100k_rr.upper) %>%
right_join(
ward_pulm_impact2 %>%
select(ward, post_effect = acf_inc100k_rr,
post_effect_lower = acf_inc100k_rr.lower,
post_effect_upper = acf_inc100k_rr.upper)
) %>%
right_join(
ward_pulm_impact3 %>%
filter(acf_period=="c. post-acf") %>%
select(ward, slope_effect = pct_change_epred_annual,
slope_effect_lower = pct_change_lower_annual,
slope_effect_upper = pct_change_upper_annual)
)
Error in select(., ward, immediate_effect = acf_inc100k_rr, immediate_effect_lower = acf_inc100k_rr.lower, :
object 'ward_impact_out' not found
Try a different way with the full distribution of posteriors
Correlation between immediate effect and post effect of ACF
Correlation between immediate effect and change in slope
Join these together with the overall estimates to make a single figure for showing impact
ward_counterf <- calcuate_counterfactual(model_data = mdata2, model=m_pulmonary_ward, population_denominator = population_without_inst_ship, grouping_var=ward)
`summarise()` has grouped output by '.draw'. You can override using the `.groups` argument.Joining with `by = join_by(year, population_without_inst_ship, .draw)`Warning: Detected an unexpected many-to-many relationship between `x` and `y`.Joining with `by = join_by(.draw)`
ward_counterf %>%
map(datatable)
$counter_post
$counter_post_overall
NA
Total pulmonary TB cases averted between 1958 and 1963
Fit the model
(Not rewritten the functions for this yet)
mdata3 <- cases_by_age_sex %>%
filter(tb_type=="Pulmonary") %>%
mutate(acf_period = case_when(year %in% c(1950:1956) ~ "a. pre-acf",
year %in% c(1957) ~ "b. acf",
year %in% c(1958:1963) ~ "c. post-acf")) %>%
mutate(year2 = year+0.5) %>%
group_by(age, sex) %>%
mutate(y_num = row_number()) %>%
ungroup()
winform_prior3 <- c(prior(normal(0, 1), class = Intercept),
#prior(gamma(0.5, 0.01), class = shape),
prior(normal(0, 1), class = b),
prior(cauchy(0,5), class="sd"),
prior(lkj(2), class="cor"))
m_age_sex <- brm(
cases ~ y_num + (acf_period)*(age*sex) + (acf_period:y_num)*(age*sex),
data = mdata3,
family = negbinomial(),
seed = 1234,
chains = 4, cores = 4,
prior = basic_prior,
backend = "cmdstanr")
Compiling Stan program...
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
In file included from /var/folders/bh/wr0g5x9j2wq46p9hm3ft_b380000gn/T/RtmpvqbQUj/model-115da1c751eb2.hpp:1:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/src/stan/model/model_header.hpp:4:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math.hpp:19:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/rev.hpp:10:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/rev/fun.hpp:198:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor.hpp:15:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor/integrate_ode_rk45.hpp:6:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:9:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint.hpp:76:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/observer_collection.hpp:23:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function.hpp:30:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function/detail/prologue.hpp:17:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/function/function_base.hpp:21:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/type_index.hpp:29:
In file included from /Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/type_index/stl_type_index.hpp:47:
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.
/
0/boost/container_hash/hash.hpp:132:33: warning: 'unary_function<const std::error_category *, unsigned long>' is deprecated [-Wdeprecated-declarations]
struct hash_base : std::unary_function<T, std::size_t> {};
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:692:18: note: in instantiation of template class 'boost::hash_detail::hash_base<const std::error_category *>' requested here
: public boost::hash_detail::hash_base<T*>
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:420:24: note: in instantiation of template class 'boost::hash<const std::error_category *>' requested here
boost::hash<T> hasher;
^
/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/boost_1.78.0/boost/container_hash/hash.hpp:551:9: note: in instantiation of function template specialization 'boost::hash_combine<const std::error_category *>' requested here
hash_combine(seed, &v.category());
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__functional/unary_function.h:23:29: note: 'unary_function<const std::error_category *, unsigned long>' has been explicitly marked deprecated here
struct _LIBCPP_TEMPLATE_VIS _LIBCPP_DEPRECATED_IN_CXX11 unary_function
^
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__config:850:41: note: expanded from macro '_LIBCPP_DEPRECATED_IN_CXX11'
# define _LIBCPP_DEPRECATED_IN_CXX11 _LIBCPP_DEPRECATED
^
-
/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/c++/v1/__config:835:49: note: expanded from macro '_LIBCPP_DEPRECATED'
# define _LIBCPP_DEPRECATED __attribute__((deprecated))
^
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
|
/
-
\
1 warning generated.
|
/
ld: warning: duplicate -rpath '/Users/petermacpherson/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
-
\
|
Start sampling
Running MCMC with 4 parallel chains...
Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 2 finished in 9.4 seconds.
Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1 finished in 9.7 seconds.
Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 finished in 9.6 seconds.
Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 finished in 9.6 seconds.
All 4 chains finished successfully.
Mean chain execution time: 9.6 seconds.
Total execution time: 10.1 seconds.
summary(m_age_sex)
Family: negbinomial
Links: mu = log; shape = identity
Formula: cases ~ y_num + (acf_period) * (age * sex) + (acf_period:y_num) * (age * sex)
Data: mdata3 (Number of observations: 224)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.47 0.08 4.31 4.62 1.00 2429 2931
y_num -0.18 0.02 -0.22 -0.13 1.00 1390 2194
acf_periodb.acf -0.01 0.25 -0.49 0.47 1.00 10433 2887
acf_periodc.postMacf 0.01 0.15 -0.30 0.30 1.00 4606 3010
age06_15 0.41 0.12 0.18 0.65 1.00 3070 2827
age16_25 1.56 0.11 1.34 1.79 1.00 2070 2939
age26_35 0.92 0.11 0.70 1.14 1.00 2537 2963
age36_45 0.24 0.12 0.02 0.47 1.00 3037 3317
age46_55 -0.34 0.13 -0.59 -0.08 1.00 2813 2851
age56_65 -0.64 0.13 -0.90 -0.38 1.00 3571 3404
age65P -0.99 0.14 -1.26 -0.71 1.00 3427 2876
sexM 0.20 0.09 0.03 0.37 1.00 2440 2716
age06_15:sexM -0.29 0.15 -0.58 0.01 1.00 3885 3215
age16_25:sexM -0.31 0.14 -0.58 -0.03 1.00 2987 3222
age26_35:sexM -0.18 0.14 -0.45 0.09 1.00 2928 2855
age36_45:sexM 0.16 0.14 -0.12 0.44 1.00 3311 3278
age46_55:sexM 0.69 0.15 0.38 0.99 1.00 2993 2969
age56_65:sexM 0.53 0.16 0.23 0.84 1.00 3329 3430
age65P:sexM 0.25 0.16 -0.06 0.56 1.00 3310 3450
acf_periodb.acf:age06_15 0.02 0.25 -0.47 0.53 1.00 8378 2757
acf_periodc.postMacf:age06_15 -0.37 0.22 -0.81 0.08 1.00 6912 2822
acf_periodb.acf:age16_25 0.03 0.25 -0.47 0.52 1.00 8484 2916
acf_periodc.postMacf:age16_25 0.23 0.20 -0.16 0.61 1.00 6331 3214
acf_periodb.acf:age26_35 0.04 0.24 -0.44 0.52 1.00 9149 2749
acf_periodc.postMacf:age26_35 0.14 0.20 -0.26 0.54 1.00 6922 3253
acf_periodb.acf:age36_45 0.05 0.25 -0.44 0.55 1.00 8346 3264
acf_periodc.postMacf:age36_45 0.10 0.20 -0.29 0.50 1.00 7240 3137
acf_periodb.acf:age46_55 0.05 0.25 -0.44 0.54 1.00 8134 2763
acf_periodc.postMacf:age46_55 0.20 0.21 -0.20 0.59 1.00 5967 3717
acf_periodb.acf:age56_65 0.04 0.24 -0.45 0.52 1.00 7812 2763
acf_periodc.postMacf:age56_65 0.05 0.21 -0.37 0.46 1.00 6009 2943
acf_periodb.acf:age65P 0.04 0.25 -0.44 0.53 1.00 8029 3377
acf_periodc.postMacf:age65P 0.11 0.22 -0.31 0.55 1.00 6647 3257
acf_periodb.acf:sexM -0.01 0.24 -0.47 0.45 1.00 9693 2911
acf_periodc.postMacf:sexM 0.17 0.17 -0.18 0.51 1.00 5883 3083
y_num:acf_periodb.acf -0.06 0.05 -0.15 0.03 1.00 2624 2474
y_num:acf_periodc.postMacf 0.00 0.03 -0.05 0.05 1.00 1687 2704
acf_periodb.acf:age06_15:sexM 0.01 0.24 -0.48 0.48 1.00 7743 2919
acf_periodc.postMacf:age06_15:sexM -0.27 0.23 -0.72 0.18 1.00 8217 3108
acf_periodb.acf:age16_25:sexM 0.01 0.24 -0.47 0.48 1.00 8818 2623
acf_periodc.postMacf:age16_25:sexM 0.12 0.23 -0.32 0.58 1.00 7353 3074
acf_periodb.acf:age26_35:sexM 0.01 0.25 -0.49 0.50 1.00 9256 2639
acf_periodc.postMacf:age26_35:sexM 0.05 0.22 -0.39 0.50 1.00 7281 2607
acf_periodb.acf:age36_45:sexM 0.00 0.25 -0.49 0.48 1.00 8703 3032
acf_periodc.postMacf:age36_45:sexM 0.01 0.22 -0.43 0.45 1.00 7249 2773
acf_periodb.acf:age46_55:sexM 0.01 0.24 -0.47 0.49 1.00 9103 2587
acf_periodc.postMacf:age46_55:sexM 0.22 0.22 -0.21 0.65 1.00 6530 3109
acf_periodb.acf:age56_65:sexM 0.02 0.25 -0.47 0.53 1.00 9465 2697
acf_periodc.postMacf:age56_65:sexM 0.10 0.22 -0.33 0.53 1.00 7499 2951
acf_periodb.acf:age65P:sexM 0.02 0.25 -0.46 0.50 1.00 9198 3071
acf_periodc.postMacf:age65P:sexM 0.12 0.23 -0.32 0.55 1.01 7371 2487
y_num:acf_perioda.preMacf:age06_15 0.06 0.03 0.00 0.13 1.00 1906 2589
y_num:acf_periodb.acf:age06_15 0.14 0.05 0.04 0.24 1.00 2804 2660
y_num:acf_periodc.postMacf:age06_15 0.08 0.02 0.03 0.12 1.00 4267 3359
y_num:acf_perioda.preMacf:age16_25 0.18 0.03 0.12 0.24 1.00 1422 2399
y_num:acf_periodb.acf:age16_25 0.25 0.05 0.15 0.34 1.00 2590 2847
y_num:acf_periodc.postMacf:age16_25 0.03 0.02 -0.02 0.07 1.00 3339 2925
y_num:acf_perioda.preMacf:age26_35 0.19 0.03 0.13 0.25 1.00 1637 2735
y_num:acf_periodb.acf:age26_35 0.30 0.05 0.21 0.40 1.00 2660 3000
y_num:acf_periodc.postMacf:age26_35 0.08 0.02 0.04 0.13 1.00 3778 3384
y_num:acf_perioda.preMacf:age36_45 0.17 0.03 0.11 0.23 1.00 1912 2636
y_num:acf_periodb.acf:age36_45 0.37 0.05 0.27 0.47 1.00 2570 2892
y_num:acf_periodc.postMacf:age36_45 0.12 0.02 0.07 0.16 1.00 3471 3230
y_num:acf_perioda.preMacf:age46_55 0.13 0.03 0.07 0.20 1.00 2070 2746
y_num:acf_periodb.acf:age46_55 0.37 0.05 0.27 0.47 1.00 2640 3015
y_num:acf_periodc.postMacf:age46_55 0.12 0.02 0.08 0.17 1.00 3329 3026
y_num:acf_perioda.preMacf:age56_65 0.09 0.04 0.01 0.16 1.00 2557 2631
y_num:acf_periodb.acf:age56_65 0.30 0.05 0.20 0.40 1.00 2875 2937
y_num:acf_periodc.postMacf:age56_65 0.12 0.02 0.07 0.17 1.00 3806 3521
y_num:acf_perioda.preMacf:age65P 0.10 0.04 0.02 0.17 1.00 2323 2883
y_num:acf_periodb.acf:age65P 0.31 0.05 0.21 0.42 1.00 2557 2529
y_num:acf_periodc.postMacf:age65P 0.12 0.02 0.08 0.17 1.00 3806 3486
y_num:acf_perioda.preMacf:sexM 0.00 0.03 -0.05 0.06 1.00 1427 2161
y_num:acf_periodb.acf:sexM -0.06 0.06 -0.17 0.05 1.00 1911 2495
y_num:acf_periodc.postMacf:sexM -0.04 0.02 -0.08 0.01 1.00 2451 2376
y_num:acf_perioda.preMacf:age06_15:sexM -0.03 0.04 -0.11 0.06 1.00 2333 2818
y_num:acf_periodb.acf:age06_15:sexM 0.04 0.06 -0.08 0.17 1.00 2246 2948
y_num:acf_periodc.postMacf:age06_15:sexM 0.06 0.03 0.00 0.12 1.00 3715 2984
y_num:acf_perioda.preMacf:age16_25:sexM -0.05 0.04 -0.13 0.02 1.00 1797 2555
y_num:acf_periodb.acf:age16_25:sexM 0.04 0.06 -0.07 0.16 1.00 2165 2572
y_num:acf_periodc.postMacf:age16_25:sexM 0.01 0.03 -0.04 0.07 1.00 2953 3168
y_num:acf_perioda.preMacf:age26_35:sexM -0.04 0.04 -0.11 0.04 1.00 1870 2546
y_num:acf_periodb.acf:age26_35:sexM 0.05 0.06 -0.07 0.17 1.00 2149 2476
y_num:acf_periodc.postMacf:age26_35:sexM 0.01 0.03 -0.04 0.07 1.00 3056 2807
y_num:acf_perioda.preMacf:age36_45:sexM 0.00 0.04 -0.07 0.08 1.00 2046 2557
y_num:acf_periodb.acf:age36_45:sexM 0.03 0.06 -0.09 0.15 1.00 2159 2970
y_num:acf_periodc.postMacf:age36_45:sexM 0.01 0.03 -0.04 0.07 1.00 3192 2893
y_num:acf_perioda.preMacf:age46_55:sexM 0.09 0.04 0.00 0.17 1.00 2176 2834
y_num:acf_periodb.acf:age46_55:sexM 0.07 0.06 -0.04 0.19 1.00 2182 2850
y_num:acf_periodc.postMacf:age46_55:sexM 0.01 0.03 -0.04 0.07 1.00 2864 2889
y_num:acf_perioda.preMacf:age56_65:sexM 0.17 0.04 0.09 0.26 1.00 2469 2962
y_num:acf_periodb.acf:age56_65:sexM 0.16 0.06 0.04 0.28 1.00 2273 3150
y_num:acf_periodc.postMacf:age56_65:sexM 0.09 0.03 0.04 0.15 1.00 3443 3481
y_num:acf_perioda.preMacf:age65P:sexM 0.15 0.05 0.07 0.24 1.00 2348 3140
y_num:acf_periodb.acf:age65P:sexM 0.19 0.06 0.07 0.31 1.00 2208 2439
y_num:acf_periodc.postMacf:age65P:sexM 0.09 0.03 0.03 0.15 1.00 2998 3054
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 98.38 33.81 52.32 180.75 1.00 1401 2192
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(m_age_sex)
pp_check(m_age_sex, type='ecdf_overlay')
Using 10 posterior draws for ppc type 'ecdf_overlay' by default.
Summarise posterior
ggsave(here("figures/s7.png"), height=10)
Saving 7 x 10 in image
nd <- mdata3 %>%
filter(year %in% c(1956:1957)) %>%
select(acf_period, y_num, age, sex)
age_sex_impact_out <-
add_epred_draws(m_age_sex,
newdata=nd) %>%
ungroup() %>%
select(acf_period, .epred, age, sex) %>%
pivot_wider(names_from = acf_period,
values_from = .epred,
values_fn = list) %>%
unnest() %>%
rename(pre_epred = 3,
post_epred = 4) %>%
mutate(acf_diff = post_epred-pre_epred,
acf_rr = post_epred/pre_epred) %>%
group_by(age, sex) %>%
mean_qi(acf_diff, acf_rr)
Warning: `cols` is now required when using `unnest()`.
ℹ Please use `cols = c(`a. pre-acf`, `b. acf`)`.
age_sex_impact_out %>%
mutate_if(is.double, ~ scales::number(x = ., accuracy = 0.01, big.mark = ",")) %>%
datatable()
`mutate_if()` ignored the following grouping variables:
f3a <- age_sex_impact_out %>%
mutate(sex = case_when(sex=="M" ~ "Male",
sex=="F" ~ "Female")) %>%
mutate(age = case_when(age=="00_05" ~ "0 to 5y",
age=="06_15" ~ "06 to 15y",
age=="16_25" ~ "16 to 25y",
age=="26_35" ~ "26 to 35y",
age=="36_45" ~ "36 to 45y",
age=="46_55" ~ "46 to 55y",
age=="56_65" ~ "56 to 65y",
age=="65+" ~ "65 & up y")) %>%
ggplot() +
geom_pointrange(aes(y=acf_rr, ymin=acf_rr.lower, ymax=acf_rr.upper, group=sex,
x=age,
colour = sex),
position = position_dodge(width = 0.25)) +
geom_hline(aes(yintercept=1), linetype=2) +
scale_colour_manual(values = c("purple", "darkorange"), name="") +
labs(x="",
y="Relative notifications (95% UI)\nACF (1957) vs. Before ACF (1956)") +
theme_ggdist() +
theme(legend.position = "bottom",
panel.border = element_rect(colour = "grey78", fill=NA))
nd <- mdata3 %>%
filter(year %in% c(1956,1958)) %>%
select(acf_period, y_num, age, sex)
#Do it with calculating incidence, then sumamrising.
age_sex_impact2 <-add_epred_draws(m_age_sex,
newdata=nd) %>%
ungroup() %>%
select(acf_period, .epred, age, sex) %>%
pivot_wider(names_from = acf_period,
values_from = .epred,
values_fn = list) %>%
unnest() %>%
rename(pre_epred = 3,
post_epred = 4) %>%
mutate(acf_diff = post_epred-pre_epred,
acf_rr = post_epred/pre_epred) %>%
group_by(age, sex) %>%
mean_qi(acf_diff, acf_rr)
Warning: `cols` is now required when using `unnest()`.
ℹ Please use `cols = c(`a. pre-acf`, `c. post-acf`)`.
age_sex_impact2 %>%
mutate_if(is.double, ~ scales::number(x = ., accuracy = 0.01, big.mark = ",")) %>%
datatable()
`mutate_if()` ignored the following grouping variables:
f3b <- age_sex_impact2 %>%
mutate(sex = case_when(sex=="M" ~ "Male",
sex=="F" ~ "Female")) %>%
mutate(age = case_when(age=="00_05" ~ "0 to 5y",
age=="06_15" ~ "06 to 15y",
age=="16_25" ~ "16 to 25y",
age=="26_35" ~ "26 to 35y",
age=="36_45" ~ "36 to 45y",
age=="46_55" ~ "46 to 55y",
age=="56_65" ~ "56 to 65y",
age=="65+" ~ "65 & up y")) %>%
ggplot() +
geom_pointrange(aes(y=acf_rr, ymin=acf_rr.lower, ymax=acf_rr.upper, group=sex,
x=age,
colour = sex),
position = position_dodge(width = 0.25)) +
geom_hline(aes(yintercept=1), linetype=2) +
scale_colour_manual(values = c("purple", "darkorange"), name="") +
labs(x="",
y="Relative notifications (95% UI)\nACF (1958) vs. Before ACF (1956)") +
theme_ggdist() +
theme(legend.position = "bottom",
panel.border = element_rect(colour = "grey78", fill=NA))
age_sex_impact3 <- mdata3 %>%
select(year, year2, y_num, acf_period, cases, age, sex) %>%
filter(year!=1957) %>%
add_epred_draws(m_age_sex) %>%
group_by(year, age, sex, acf_period) %>%
mean_qi(.epred) %>%
ungroup() %>%
mutate(n_years = length(year), .by=acf_period) %>%
summarise(pct_change_epred_overall = (((last(.epred) - first(.epred))/first(.epred))),
pct_change_lower_overall = (((last(.lower) - first(.lower))/first(.lower))),
pct_change_upper_overall = (((last(.upper) - first(.upper))/first(.upper))),
pct_change_epred_annual = (((last(.epred) - first(.epred))/first(.epred))/n_years),
pct_change_lower_annual = (((last(.lower) - first(.lower))/first(.lower))/n_years),
pct_change_upper_annual = (((last(.upper) - first(.upper))/first(.upper))/n_years),
.by = c(acf_period, age, sex)) %>%
distinct()
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.
age_sex_impact3 %>%
mutate_if(is.double, percent) %>%
datatable()
f3c <- age_sex_impact3 %>%
mutate(sex = case_when(sex=="M" ~ "Male",
sex=="F" ~ "Female")) %>%
mutate(age = case_when(age=="00_05" ~ "0 to 5y",
age=="06_15" ~ "06 to 15y",
age=="16_25" ~ "16 to 25y",
age=="26_35" ~ "26 to 35y",
age=="36_45" ~ "36 to 45y",
age=="46_55" ~ "46 to 55y",
age=="56_65" ~ "56 to 65y",
age=="65+" ~ "65 & up y")) %>%
ggplot() +
geom_pointrange(aes(y=pct_change_epred_annual, ymin=pct_change_lower_annual, ymax=pct_change_upper_annual, group=acf_period,
x=age,
colour = acf_period), size=0.1) +
scale_y_continuous(labels =percent) +
facet_grid(.~sex) +
coord_flip() +
scale_colour_manual(values = c("#DE0D92", "#4D6CFA")) +
labs(x="",
y="Mean annual rate of change in case notification rate (95% UI)\n Before ACF (1950-1956) vs. after ACF (1958-1963)",
colour="") +
theme_ggdist() +
theme(panel.border = element_rect(colour = "grey78", fill=NA))
f3c
Join together for Figure 2.
(f3a / f3b / f3c) + plot_annotation(tag_levels = "A")
ggsave(here("figures/f3.png"), height=12)
Saving 7.29 x 12 in image
counter_post_overall_age_sex <-
left_join(counterfact_overall_age_sex, post_change_overall_age_sex) %>%
mutate(cases_averted = .epred_counterf-.epred,
pct_change = (.epred - .epred_counterf)/.epred_counterf) %>%
group_by(age, sex) %>%
mean_qi(cases_averted, pct_change) %>%
ungroup() %>%
mutate(year = "Overall (1958-1963)")
Joining with `by = join_by(age, sex, .draw)`
(Very much a work in progress!)
plot_counterfactual(model=m_pulmonary_division_prior, model_data=mdata3, population_denominator = population_without_inst_ship, outcome = cases, grouping_var = division, division) + scale_y_log10()
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Make a table of counterfactual effects for the manuscript